
############################################
############################################
# previous: 5-VAP_CurrentPlotsThresh########
############################################

###### Summary ############################################################################################
#This script unzips al future layers into the future layer directory and writes shell files that project habitat suitability for each 
#modelling unit to al future climate scenarios.
###########################################################################################################

############################################
# next: 7-VAP_FutSum_Thresh_Plot############
############################################

#######################################################################################################
##### MaxEnt Future Projections after variable selection ##############################################
#######################################################################################################

#load packages:
library(R.utils)
library(raster)
library(maptools)
library(rgdal)
library(stringr)
library(sf)


#######################################################################################################
#### the below settings may change depending on the project ###########################################
#######################################################################################################

#set project name:
projectName<-"WHOSnakes"

#define selection criteria of interest for further analyses:
crits<-c("permutation.importance")   #c("permutation.importance","contribution","Test.gain.without","Test.gain.with.only","AUC.without","AUC.with.only")
critThresh<-c(1)                     #c(1,1,1,10,1,75)

#define layer directory to use for data extractions and/or projections:


#######################################################################################################
##### the below settings should all stay the same in most cases #######################################
#######################################################################################################

#define dirs:
MyDir<-paste("/home/jc217070/Maxent",sep="")
maxent.jar <- paste(MyDir,"/maxent.jar",sep="")

LayersDir<-paste(MyDir, "/Layers",sep="")
FutLayersDir<-paste("/home/jc217070/GIS_and_raw_Layers/final_layers/WorldClim2.1/future")
scenarios<-list.dirs(path=paste(FutLayersDir),full.names=FALSE,recursive=FALSE)

OutDirFinal<-paste(MyDir,"/Outputs3_Final",sep="")

PbsDir <- paste(MyDir, "/tmp.pbs", sep="")
setwd(PbsDir)

#read species reference list:
DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
spp<-DatSum$gen_sp_subtax
mapspp<-unique(DatSum$ModUnit)
#mapsppB<-unique(DatSum$ModUnit[DatSum$Priority==0 | DatSum$Priority==1])
#mapsppB<-unique(DatSum$ModUnit[DatSum$Redo==1])
mapsppB<-mapspp
#mapsppB<- mapspp[!(mapspp %in% (c("Echis ocellatus", "Echis leucogaster","Naja senegalensis","Naja savannula","Bitis arietans","Dendroaspis polylepis")))]
#mapsppB<- (c("Echis ocellatus", "Echis leucogaster","Naja senegalensis","Naja savannula","Bitis arietans","Dendroaspis polylepis"))
#mapsppB<-unique(DatSum$ModUnit[DatSum$Subclade %in% c("Echis","Naja sensu stricto","Bitis","Dendroaspis; Ophiophagus")])




#make tet file with qsub commands list for all species' shell files
qsub.list.file <- '/home/jc217070/Maxent/R_scripts/qsubs-6_VAP.txt'

#################################################################
# set up layer dir for future projections #######################
#################################################################

C_WCfiles<-list.files(path=LayersDir,pattern='wc2.1_')
C_Allfiles<-list.files(path=LayersDir,pattern='.asc')
C_nonWCfiles<-C_Allfiles[!(C_Allfiles %in% C_WCfiles)]

#for (scen in scenarios[c(1,3,5,7,9,11,13,15,17,19,21,23,25,27,29)]){
# for (scen in c("INM-CM5-0_ssp585_2041-2060","MPI-ESM1-2-LR_ssp585_2041-2060","CanESM5-CanOE_ssp585_2041-2060",
#                "EC-Earth3-Veg-LR_ssp585_2041-2060","FIO-ESM-2-0_ssp585_2041-2060","CMCC-ESM2_ssp585_2041-2060",
#                "INM-CM4-8_ssp585_2041-2060")){
for (scen in c("INM-CM5-0_ssp585_2041-2060","MPI-ESM1-2-LR_ssp585_2041-2060","CanESM5-CanOE_ssp585_2041-2060",
                 "EC-Earth3-Veg-LR_ssp585_2041-2060","FIO-ESM-2-0_ssp585_2041-2060","CMCC-ESM2_ssp585_2041-2060",
                 "INM-CM4-8_ssp585_2041-2060","INM-CM5-0_ssp585_2081-2100","MPI-ESM1-2-LR_ssp585_2081-2100","CanESM5-CanOE_ssp585_2081-2100",
                 "EC-Earth3-Veg-LR_ssp585_2081-2100","FIO-ESM-2-0_ssp585_2081-2100","CMCC-ESM2_ssp585_2081-2100",
                 "INM-CM4-8_ssp585_2081-2100")){
    
 for (file in C_nonWCfiles){
   if(!file.exists(paste(FutLayersDir,"/",scen,"/",file,sep=""))){
   #copy non WC files to new location
   print(paste(Sys.time(), "copying", file, "for scenario",scen))
   file.copy(from=paste(LayersDir,"/",file,sep=""), 
             to=paste(FutLayersDir,"/",scen,"/",file,sep=""),
             overwrite=TRUE,recursive=FALSE,
             copy.mode = TRUE, copy.date = FALSE)
 }}
  
  F_WCfiles<-list.files(path=paste(FutLayersDir,"/",scen, sep=""),pattern='.asc.gz')
  F_WCfiles<-unique(gsub(".gz","",F_WCfiles))
  
  for (file in F_WCfiles){
  #unzip future WC files
  old_name<-paste(gsub("wc2.1_30s_bioc_",paste("wc2.1_"),gsub(scen,"",gsub("_No","bio",gsub("_RAD_","RAD",gsub("_RH_","RH",gsub("of19","",gsub(".asc.gz",".asc",file))))))),sep="")
  if(old_name %in% C_WCfiles & !file.exists(paste(FutLayersDir,"/",scen,"/",old_name,sep=""))){
    print(paste(Sys.time(), "unzipping", file, "for scenario",scen))
    gunzip(filename=paste(FutLayersDir,"/",scen,"/",file,".gz",sep=""),
         destname=paste(FutLayersDir,"/",scen,"/",old_name,sep=""),
         remove=FALSE, overwrite=TRUE)
  }}
  #}
  
  
##############################################################

#for (scen in scenarios){
  
  for (sp in mapspp){

  ##########################
  if (sp %in% mapsppB){
  ##########################
    
  n<-which(DatSum$ModUnit==sp) #gets row number for current species
  sp<-gsub(" ","_",sp)
  
  for (crit in crits){ #really only permutation importance
    
      ###### shellfile future #######
      # create the shell file for each species model job & submit to HPC
      
      EnvLayDir<-file.path(FutLayersDir,scen)
      SpOutDirFin<-paste(OutDirFinal, "/", crit,"/ModUnit_",sp, sep="")
      
      if(gsub("2081-2100","",scen)==scen){
        year<-"2050"
        yearlong<-"2041-2060"
        
      }else{
        year<-"2090"
        yearlong<-"2081-2100"
      }
      
      dir.create(paste(SpOutDirFin, "/future", sep=""))
      outfilename <-   paste(SpOutDirFin, "/future/ModUnit_",sp,"_",scen, sep="")
      
      
      
      allfilenamesASC <-   paste(SpOutDirFin, "/future/ModUnit_",sp,"_",unique(gsub("_2081-2100","",gsub("_2041-2060","",scenarios))),"_",yearlong,".asc", sep="")
      allfilenamesGZ <-   paste(SpOutDirFin, "/future/ModUnit_",sp,"_",unique(gsub("_2081-2100","",gsub("_2041-2060","",scenarios))),"_",yearlong,".asc.gz", sep="")
      

      filesizesASC<-3500000000
      for(file in allfilenamesASC){
        if(file.exists(file)){
        filesize<-as.numeric(file.info(file)[1])
        filesizesASC<-c(filesizesASC,filesize)
        }
      }

      filesizesGZ<-65000000
      for(file in allfilenamesGZ){
        if(file.exists(file)){
          filesize<-as.numeric(file.info(file)[1])
          filesizesGZ<-c(filesizesGZ,filesize)
        }
      }
      
      
      
      
      
      
      if((!file.exists(paste(outfilename,".asc",sep="")) & !file.exists(paste(outfilename,".asc.gz",sep=""))) |
         (!file.exists(paste(outfilename,"_limitFactor.asc",sep="")) & !file.exists(paste(outfilename,"_limitFactor.asc.gz",sep="")))|

         (file.exists(paste(outfilename,".asc",sep="")) & !file.exists(paste(outfilename,".asc.gz",sep="")) & as.numeric(file.info(paste(outfilename,".asc",sep=""))[1]) < (max(filesizesASC)-300000000))|
         (file.exists(paste(outfilename,".asc.gz",sep="")) &!file.exists(paste(outfilename,".asc",sep=""))& as.numeric(file.info(paste(outfilename,".asc.gz",sep=""))[1])< (max(filesizesGZ)- 45000000))|

          (file.exists(paste(outfilename,"_limitFactor.asc",sep="")) & !file.exists(paste(outfilename,"_limitFactor.asc.gz",sep="")) & as.numeric(file.info(paste(outfilename,"_limitFactor.asc",sep=""))[1])  <3000000000)|
          (file.exists(paste(outfilename,"_limitFactor.asc.gz",sep="")) & !file.exists(paste(outfilename,"_limitFactor.asc",sep="")) & as.numeric(file.info(paste(outfilename,"_limitFactor.asc.gz",sep="")))[1]< 7000000)){

      
      
      
      if (file.exists(paste(OutDirFinal,"/",crit,"/ModUnit_", sp, "/", sp, ".html", sep=""))){ #safety, make sure again that output exists
        
        Sys.sleep(0.1)
        print(paste(Sys.time(), "projecting", sp, crit, "model to future scenario", scen))
        
        shell.file.name <- paste(PbsDir, "/ModUnit_", sp,"_",crit, "_",scen,"_proj.sh", sep="")
        
        shell.file <- file(shell.file.name, "w")
        
        #complusory bits for PBS
        cat('#!/bin/bash\n', file=shell.file)
        cat('#PBS -j oe\n', file=shell.file)                         # combine stdout and stderr into one file
        cat('#PBS -m ae\n', file=shell.file)                         # combine stdout and stderr into one file
        cat('#PBS -l walltime=10:00:00\n', sep='', file=shell.file)          # define walltime hh:mm:ss !!!!!!!!!!!!!!!!!!
        cat('#PBS -l select=1:ncpus=1:mem=70gb\n', sep='', file=shell.file)                  # allocate memory !!!!!!!!!!!!!!!!!
        
        cat('cd $PBS_O_WORKDIR\n', file=shell.file)
        cat('shopt -s expand_aliases\n', file=shell.file)
        cat('source /etc/profile.d/modules.sh\n', file=shell.file)            # need to access modules for java
        cat('echo "Job identifier is $PBS_JOBID"\n', file=shell.file)         # need this to run R
        cat('echo "Working directory is $PBS_O_WORKDIR"\n', file=shell.file)  # need this to run R
        cat('module load java\n', file=shell.file)                            # need to load java for maxent
        
        # project maxent model
        outfilename <-   paste(SpOutDirFin, "/future/ModUnit_",sp,"_",scen, sep="")
        cat('java -cp ',maxent.jar, ' density.Project ',SpOutDirFin, '/', sp, '.lambdas ',EnvLayDir, ' ', outfilename, '.asc nowriteclampgrid noextrapolate dontextrapolate \n', sep="", file=shell.file)

        #create limiting factor map:
        limitfilename <- paste(SpOutDirFin, "/future/ModUnit_",sp,"_",scen,"_limitFactor", sep="")
        cat('java -cp ',maxent.jar,' density.tools.LimitingFactor ' , SpOutDirFin, '/', sp, '.lambdas ', LayersDir, ' ', limitfilename, '.asc', sep="",file=shell.file)

        # zip output ascis:
        #cat('gzip -f ', outfilename, '.asc\n', sep="", file=shell.file)
        #cat('gzip -f ', limitfilename, '.asc\n', sep="", file=shell.file)

        close(shell.file)
        
        write(paste("qsub ", shell.file.name, sep=""), file=qsub.list.file, append=T) 
        
        #system(paste("qsub -p -100 ", shell.file.name, sep="")) #set low priority so that I can still run other jobs while all these are running
        #Sys.sleep(5) # wait 5 sec between job submissions
     }
        } #else {
    # if (file.exists(paste(outfilename,".asc",sep=""))){
    #   Sys.sleep(0.1)
    #   print(paste(Sys.time(), "zipping", sp, scen, "future projection"))
    #   gzip(filename=paste(outfilename,".asc",sep=""),overwrite=TRUE, remove=TRUE)
    #   if (file.exists(paste(outfilename,"_limitFactor.asc",sep=""))){
    #     print(paste(Sys.time(), "zipping", sp, scen, "limit factor raster"))
    #   gzip(filename=paste(outfilename,"_limitFactor.asc",sep=""),overwrite=TRUE, remove=TRUE)
    # }}}
  }
}
}
}

#cd /home/jc217070/Maxent/tmp.pbs
#bash /home/jc217070/Maxent/R_scripts/qsubs-6_VAP.txt

############################################
############################################
# next: 7-VAP_FutSum_Thresh_Plot############
############################################




